#*********************************************************#
#                                                         #
#                                                         #
#                 REPLICATION FILES FOR:                  #
#                                                         #
#                 Barcel� and Labzina                     #
#                                                         #
# Do Islamic State's Deadly Attacks Disengage, Deter, or  #
#                 Mobilize Supporters?                    #
#                                                         #
#         British Journal of Political Science            #
#                                                         #
#                                                         #         
#*********************************************************#


#***************#
# Load Packages #
#***************#

library(foreign) 
library(boot)
library(stargazer)
library(MASS)
library(Hmisc)
library(effects)
library(xtable) 
library(lme4)
library(gridExtra)
library(ggplot2);
library(car)
library(plm)
library(reshape)
library(lmtest)
library(rdd)
library(rddtools)

#***************#
#   Load Data   #
#***************#

#Set Working Directory

path <- ""

setwd(path)

re_data <- read.csv2("re_data.csv", sep = ",")
fe_data <- read.csv2("fe_data.csv", sep = ",")

re_data[,6:16] <- sapply(sapply(re_data[,6:16],as.character), as.numeric)
fe_data[,6:16] <- sapply(sapply(fe_data[,6:16],as.character), as.numeric)

###############
####TABLE 1####
###############

rd_brussels_1d <- RDestimate(followers_count ~ time_brussels, cutpoint= -0.5, bw = c(1,1,1), data= fe_data[fe_data$time_brussels > -1.1 & fe_data$time_brussels < 0.1,]) #1 day

round(rd_brussels_1d$est, 1)
round(rd_brussels_1d$se, 1)
round(rd_brussels_1d$p/2, 2)
rd_brussels_1d$obs

rd_brussels <- RDestimate(followers_count ~ time_brussels, cutpoint= -0.5, data=fe_data)

rd_brussels$bw
round(rd_brussels$est, 1)
round(rd_brussels$se, 1)
round(rd_brussels$p/2,2)
rd_brussels$obs

rd_nice_1d <- RDestimate(followers_count ~ time_nice, cutpoint= -0.5, bw = c(1,1,1), data= fe_data[fe_data$time_nice > -1.1 & fe_data$time_nice < 0.1,]) #1 day

round(rd_nice_1d$est, 2)
round(rd_nice_1d$se, 2)
round(rd_nice_1d$p/2, 2)
rd_nice_1d$obs

rd_nice <- RDestimate(followers_count ~ time_nice, cutpoint= -0.5, data=fe_data)

rd_nice$bw
round(rd_nice$est, 2)
round(rd_nice$se, 2)
round(rd_nice$p/2, 2)
rd_nice$obs

###############
####TABLE 2####
###############

#PANEL A

#Column 1.1
summary(tab2_col1.1 <- lmer(lnfollowers_count ~ Deaths_EurUS50_00 + Deaths_nonEurUS50_00 + (1 | time) + (1 | account_id), data=re_data))

#Column 1.2
summary(tab2_col1.2 <- plm(lnfollowers_count ~ Deaths_EurUS50_00 + Deaths_nonEurUS50_00, data= fe_data, model="within", index=c("account_id", "time")))

coeftest(tab2_col1.2, vcov=vcovHC(tab2_col1.2,type="HC0",cluster="group"))

#Column 1.3
summary(tab2_col1.3 <- lmer(log(followers_count+1) ~ Deaths_EurUS50_00 + Deaths_nonEurUS50_00 + Nreports_0000 + nsuspensions_000 + (1 | time) + (1 | account_id), data=re_data))

#Column 1.4
summary(tab2_col1.4 <- plm(lnfollowers_count ~ Deaths_EurUS50_00 + Deaths_nonEurUS50_00 + Nreports_0000 + nsuspensions_000, data= fe_data, model="within", index=c("account_id", "time")))

coeftest(tab2_col1.4, vcov=vcovHC(tab2_col1.4,type="HC0",cluster="group"))

#PANEL B

#Column 2.1
summary(tab2_col2.1 <- lmer(lnfollowers_count ~ Deaths_EurUS75_00 + Deaths_nonEurUS75_00 + (1 | time) + (1 | account_id), data=re_data))

#Column 2.2
summary(tab2_col2.2 <- plm(lnfollowers_count ~ Deaths_EurUS75_00 + Deaths_nonEurUS75_00, data= fe_data, model="within", index=c("account_id", "time")))

coeftest(tab2_col2.2, vcov=vcovHC(tab2_col2.2,type="HC0",cluster="group"))

#Column 2.3
summary(tab2_col2.3 <- lmer(lnfollowers_count ~ Deaths_EurUS75_00 + Deaths_nonEurUS75_00 + Nreports_0000 + nsuspensions_000 + (1 | time) + (1 | account_id), data=re_data))

#Column 2.4
summary(tab2_col2.4 <- plm(lnfollowers_count ~ Deaths_EurUS75_00 + Deaths_nonEurUS75_00 + Nreports_0000 + nsuspensions_000, data= fe_data, model="within", index=c("account_id", "time")))

coeftest(tab2_col2.4, vcov=vcovHC(tab2_col2.4,type="HC0",cluster="group"))

#PANEL C

#Column 3.1
summary(tab2_col3.1 <- lmer(lnfollowers_count ~ Deaths_EurUS100_00 + Deaths_nonEurUS100_00 + (1 | time) + (1 | account_id), data=re_data))

#Column 3.2
summary(tab2_col3.2 <- plm(lnfollowers_count ~ Deaths_EurUS100_00 + Deaths_nonEurUS100_00, data= fe_data, model="within", index=c("account_id", "time")))

coeftest(tab2_col3.2, vcov=vcovHC(tab2_col3.2,type="HC0",cluster="group"))

#Column 3.3
summary(tab2_col3.3 <- lmer(lnfollowers_count ~ Deaths_EurUS100_00 + Deaths_nonEurUS100_00 + Nreports_0000 + nsuspensions_000 + (1 | time) + (1 | account_id), data=re_data))

#Column 3.4
summary(tab2_col3.4 <- plm(lnfollowers_count ~ Deaths_EurUS100_00 + Deaths_nonEurUS100_00 + Nreports_0000 + nsuspensions_000, data= fe_data, model="within", index=c("account_id", "time")))

coeftest(tab2_col3.4, vcov=vcovHC(tab2_col3.4,type="HC0",cluster="group"))

###############
####TABLE 3####
###############

#Simulated values from models#

###############
####TABLE 4####
###############

re_data_cinc <- read.csv2("re_data_cinc.csv") 
fe_data_cinc <- read.csv2("fe_data_cinc.csv")

summary(tab4_col1 <- lm(lnfollowers_count ~ scores*(Deaths_EurUS50_00 + Deaths_nonEurUS50_00), data=re_data_cinc))

summary(tab4_col2 <- plm(lnfollowers_count ~ scores*(Deaths_EurUS50_00 + Deaths_nonEurUS50_00), data= fe_data_cinc, model="within", index=c("X.account_id", "time")))

coeftest(tab4_col2, vcov=vcovHC(tab4_col2,type="HC0",cluster="group"))

summary(tab4_col3 <- lm(lnfollowers_count ~ scores*(Deaths_EurUS50_00 + Deaths_nonEurUS50_00) + Nreports_000 + nsuspensions_00, data=re_data_cinc))

summary(tab4_col4 <- plm(lnfollowers_count ~ scores*(Deaths_EurUS50_00 + Deaths_nonEurUS50_00) + Nreports_000 + nsuspensions_00, data= fe_data_cinc, model="within", index=c("X.account_id", "time")))
coeftest(tab4_col4, vcov=vcovHC(tab4_col4,type="HC0",cluster="group"))

###############
####FIGURE 1###
###############

attacks <- read.csv2("attacks.csv", sep = ',')

df <- data.frame(date=as.Date(attacks$�..date, "%d.%m.%Y"),v1=as.numeric(as.character(attacks$Deaths_nonEurUS50)), v2=as.numeric(as.character(attacks$Deaths_EurUS50)))

newdf <- melt(df,'date')

p50 <- ggplot(newdf,aes(x=date,y=value, fill=variable)) + theme_bw() +
  geom_bar(stat="identity") +
  theme(legend.title=element_blank()) +
  xlab("Date") + ylab("Number of deaths (50% discount rate)") +
  scale_fill_manual(values=c("v1" = 'gray75', 'v2' = "black"), labels =c('World', 'US and Europe')) +
  theme(legend.position=c(.10, .90)) + theme(text = element_text(size=20))
p50

###############
####FIGURE 2###
###############

par(pty="m", oma=c(1, 2, 0,0))
plot(RDestimate(followers_count ~ time_brussels, cutpoint= -0.5, data=fe_data[fe_data$time_brussels < 15,]))
title(ylab="Number of followers in IS-related Accounts", xlab="Time to Attack (in days)")

par(pty="m", oma=c(1, 2, 0,0))
plot(RDestimate(followers_count ~ time_nice, cutpoint= -0.5, data=fe_data[fe_data$time_nice > -15,]))
title(ylab="Number of followers in IS-related Accounts", xlab="Time to Attack (in days)")


###############
####FIGURE 3###
###############

eff1 <- allEffects(tab2_col1.1, confidence.level=.95, xlevels=list(Deaths_EurUS50_00=c(seq(0,1,0.01))), offset = 0, transformation=list(link=log, inverse=exp))

Europe <- as.data.frame(rep(1,101))
predframe1 <- cbind(eff1$Deaths_EurUS50_00$fit, eff1$Deaths_EurUS50_00$x,
                    eff1$Deaths_EurUS50_00$se, eff1$Deaths_EurUS50_00$lower,
                    eff1$Deaths_EurUS50_00$upper, Europe)

Europe <- as.data.frame(rep(0,101))
eff2 <- allEffects(tab2_col1.1, confidence.level=.95, xlevels=list(Deaths_nonEurUS50_00=c(seq(0,1,0.01))),
                   transformation=list(link=log, inverse=exp))

predframe2 <- cbind(eff2$Deaths_nonEurUS50_00$fit, eff2$Deaths_nonEurUS50_00$x,
                    eff2$Deaths_nonEurUS50_00$se, eff2$Deaths_nonEurUS50_00$lower,
                    eff2$Deaths_nonEurUS50_00$upper, Europe)
colnames(predframe1) <- c("fit", "x", "se", "LL", "UL", "Eur")
colnames(predframe2) <- c("fit", "x", "se", "LL", "UL", "Eur")

predframe <- rbind(predframe1, predframe2)

predframe$expfit <- exp(predframe$fit)
predframe$expLL <- exp(predframe$LL)
predframe$expUL <- exp(predframe$UL)
predframe$Eur <- as.factor(predframe$Eur)

p <- ggplot(predframe, aes(x = x, y = expfit, group=Eur)) +
  theme_bw() +
  geom_point() +
  geom_line() +
  geom_ribbon(aes(ymax = expUL, ymin = expLL, fill=Eur),
              alpha = 0.7, color=NA) + 
  xlab("Cumulative Number of Victims in IS Terrorist Attacks") +
  ylab("Expected Number of Followers in IS-related Accounts") +
  theme(legend.position=c(.10, .08)) + theme(text = element_text(size=20)) + 
  scale_fill_manual(values=c("gray75", "black"), labels = c("Elsewhere", "US and Europe")) +
  guides(fill=guide_legend(title="Attack Location"))
p
